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Abstract 

This  paper  addresses  the  development  of  inverse  compensation  techniques  for  a  class  of 
ferromagnetic  transducers  including  magnetostrictive  actuators.  If  unaccommodated,  the  hys¬ 
teresis  and  nonlinear  dynamics  can  produce  severe  loss  of  control  authority  and  potential 
instabilities  when  the  actuators  are  incorporated  in  control  design.  In  this  work,  hysteresis 
is  modeled  through  the  domain  wall  theory  originally  proposed  by  Jiles  and  Atherton.  This 
model  is  based  on  the  quantification  of  the  energy  required  to  translate  domain  walls  pinned 
at  inclusions  in  the  material  with  the  magnetization  at  a  given  field  level  specified  through 
the  solution  of  an  ordinary  differential  equation.  A  complementary  differential  equation  is 
then  employed  to  compute  the  inverse  which  can  be  used  to  compensate  for  hysteresis  and 
nonlinear  dynamics  in  control  design.  The  performance  of  the  inverse  compensator  and  its 
employment  in  LQR  control  design  are  illustrated  through  numerical  examples. 
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1  Introduction 


Increased  demands  on  control  transducers  in  combination  with  novel  material  designs  have 
led  to  the  development  of  a  class  of  ferromagnetic  transducers,  including  those  based  on  mag- 
netostrictive  materials,  which  deliver  high  level  strains  and  forces  but  often  exhibit  significant 
hysteresis  and  nonlinear  dynamics.  While  magnetostrictive  transducers  can  be  utilized  in  lin¬ 
ear  regimes  by  maintaining  low  input  levels,  control  requirements  often  dictate  that  they  be 
driven  at  the  high  drive  levels  where  hysteresis  and  nonlinearities  are  inherent  to  the  actua¬ 
tor  dynamics.  Furthermore,  the  advantages  of  these  transducers  over  traditional  actuators  is 
typically  realized  at  high  drive  levels.  Hence  it  is  necessary  to  design  control  methods  which 
accommodate  the  observed  actuator  hysteresis  and  nonlinearities. 

While  numerous  mechanisms  can  produce  hysteresis  and  nonlinearities  in  magnetostrictive 
actuators,  a  primary  source  is  attributed  to  ferromagnetic  domain  mechanisms.  This  produces 
sigmoid  hysteresis  curves  of  the  type  illustrated  in  Figure  1  for  the  magnetostrictive  material 
Terfenol-D  and  is  the  form  of  hysteresis  that  we  consider  here. 

To  illustrate  the  issues  pertinent  to  control  design  using  high  performance  actuators,  we 
consider  a  plant  with  nonlinear  input  u  provided  by  the  actuator  as  indicated  in  Figure  2.  The 
input  to  the  actuator  is  denoted  by  v  while  the  dynamics  of  the  actuator  are  represented  by 
iV(-).  We  note  that  the  map  N(-)  is  typically  both  nonlinear  and  a  multivalued  function  of  the 
input  level  due  to  the  hysteresis.  It  is,  however,  a  well-defined  function  of  v.  Finally,  we  let  up 
denote  the  input  specified  to  obtain  the  control  objective  in  the  absence  of  the  nonlinearities 
and  hysteresis  N(-). 

For  systems  exhibiting  even  mild  hysteresis,  linear  control  methods  are  typically  ineffective 
due  to  incurred  phase  shifts  and  unmodeled  energy  loss  in  the  hysteresis  loop.  To  address  this, 
we  consider  the  construction  of  exact  and  approximate  inverses  for  the  actuator  nonlinearities 
and  hysteresis.  This  entails  the  computation  of  a  map  iV_1(-)  or  approximate  map  iV_1(-) 
such  that  up  =  IV-1  (v)  ~  Ar_1(n).  For  linear  control  implementation,  up  is  filtered  by  iV_1(-) 
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Figure  1.  Relationship  between  the  magnetic  field  H  and  magnetization  M  for  Terfenol-D 
(Experimental  data  collected  at  Iowa  State  University). 
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Figure  2.  (a)  Plant  with  input  actuator  hysteresis  TV(-);  (b)  Inverse  compensator  jV_1(-)  for 
the  hysteresis. 

before  input  to  the  filter.  In  this  manner,  one  compensates  for  the  nonlinearity  and  the  final 
signal  u  —  N(N^1(up))  is  that  prescribed  by  the  control  law.  This  process  is  illustrated  in 
Figure  2. 

We  emphasize  that  the  process  of  inverse  compensation  for  nonlinearities  is  well  known 
and  the  reader  is  referred  to  [1]  for  a  comprehensive  discussion  in  the  context  of  general  actu¬ 
ator  nonlinearities  and  hysteresis.  The  models  employed  in  [1]  to  characterize  the  hysteresis 
N(-)  and  inverse  jV_1(-)  are  based  on  piecewise  linear  Preisach  approximations  and  hence 
are  phenomenological  in  nature.  The  contribution  of  this  paper  lies  in  the  construction  of 
a  hysteresis  model  and  inverse  for  magnetostrictive  actuators  which  is  based  upon  energy 
principles.  This  yields  a  compensation  filter  which  is  easily  constructed  and  can  be  readily 
updated  to  accommodate  changing  operating  conditions. 

The  model  for  the  actuator  hysteresis  process  N(-)  is  based  upon  the  theory  of  Jiles  and 
Atherton  [2,  3].  This  theory  is  based  upon  the  tenant  that  ferromagnetic  hysteresis  is  primarily 
due  to  the  energy  required  to  bend  and  translate  domain  walls  pinned  at  inclusions  in  the 
material.  By  formulating  the  energy  of  domain  wall  characteristics  in  terms  of  bulk  material 
properties,  Jiles  and  Atherton  obtained  a  model  which  requires  only  five  parameters.  The 
construction  of  the  model  is  further  enhanced  by  the  physical  nature  of  the  parameters  and 
the  fact  that  certain  parameters  can  be  directly  specified  from  measured  data.  This  model 
has  subsequently  been  extended  to  magnetostrictive  transducers  [4,  5]  and  it  is  in  this  setting 
that  we  consider  the  development  of  an  inverse  jV_1(-).  We  note  that  while  the  model  and 
inverse  compensator  are  illustrated  in  the  context  of  LQR  feedback  control,  the  techniques 
are  equally  applicable  when  designing  feedforward  or  adaptive  control  systems.  Again,  the 
reader  is  referred  to  [1]  for  a  comprehensive  discussion  of  adaptive  control  design  for  systems 
exhibiting  input  hysteresis  once  an  exact  or  approximate  inverse  has  been  determined. 

Section  2  provides  a  summary  of  the  hysteresis  model  and  the  inverse  model  is  developed 
in  Section  3.  The  use  of  the  inverse  compensation  scheme  for  optimal  control  design  in  a 
structural  system  is  illustrated  in  Section  4.  For  that  development,  we  consider  magnetostric¬ 
tive  transducers  employed  as  control  actuators  for  attenuating  vibrations  in  a  thin  beam.  The 
hysteretic  actuator  inputs  N(-)  are  incorporated  in  a  thin  beam  model  and  the  dynamics  are 
approximated  to  obtain  a  finite  dimensional  control  system.  Finally,  numerical  examples  are 
provided  to  illustrate  the  performance  of  an  LQR  control  with  the  inverse  iV_1(-)  used  to 
compensate  for  the  actuator  hysteresis  and  nonlinearities. 
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2  Hysteresis  Model 

The  model  described  here  is  based  on  the  tenant  that  hysteresis  in  ferromagnetic  materials 
is  due  to  the  impeded  movement  of  domain  walls  which  are  pinned  at  inclusions  or  defects 
in  the  material.  In  such  materials,  changes  in  magnetization  are  attributed  to  the  nucleation 
and  growth  of  domains  through  domain  wall  movement  (note  that  domains  denote  regions 
in  which  moments  are  highly  aligned  while  domain  walls  are  the  transition  regions  between 
domains).  For  materials  which  are  free  from  inclusions,  domain  wall  movement  is  unimpeded 
and  the  magnetization  is  anhysteretic  (hysteresis- free).  This  situation  is  idealized,  however, 
since  defects  are  unavoidable  (e.g.,  carbons  in  steel)  and  in  many  cases,  incorporated  in  the 
material  to  achieve  the  desired  stoichiometry  (e.g.,  Dysprosium  in  Terfenol-D).  These  defects 
or  inclusions  provide  pinning  sites  for  domain  walls  due  to  the  reduction  in  energy  which  occurs 
when  the  domain  intersects  the  pinning  site.  For  low  magnetic  field  variations  about  some 
equilibrium  value,  the  walls  remain  pinned  and  the  magnetization  is  reversible.  This  motion 
becomes  irreversible  at  higher  field  levels  due  to  wall  intersections  with  remote  pinning  sites. 
The  model  quantifies  hysteresis  through  the  characterization  of  the  anhysteretic  magnetization 
Man ,  the  reversible  magnetization  Mrev  and  the  irreversible  magnetization  Mirr. 

The  anhysteretic  magnetization  at  a  point  in  the  material  is  dependent  on  the  effective 
magnetic  field,  the  saturation  magnetization  and  the  thermal  energy  of  the  sample.  Under 
the  assumption  of  constant  stress  <r0,  the  effective  field  can  be  expressed  as 


Hef  f  =  H  +  aM 


where  H  denotes  the  external  field  input  to  the  transducer,  M  denotes  the  magnetization  and 
the  parameter  a  quantifies  the  effects  of  interdomain  coupling  and  magnetoelastic  domain 
interactions.  Two  commonly  employed  models  for  the  anhysteretic  are  the  Langevin  expression 


Man  =  Ms 


coth 


(1) 


and  the  Ising  spin  relation 

Man  =  Ms tanh  •  (2) 

Both  models  are  derived  using  arguments  from  Boltzmann  statistics  with  the  differences  de¬ 
pendent  upon  the  assumptions  concerning  possible  domain  orientations.  Taylor  expansions  of 
the  hyperbolic  cotangent  and  tangent  functions  reveals  that  the  models  agree  through  third 
order  terms  with  the  Ising  spin  model  producing  steeper  slopes  due  to  the  magnitude  of 
the  higher  order  terms.  The  constant  a  is  given  by  a  =  ■  where  ks  is  Boltzmann’s  con¬ 

stant,  ksT  represents  the  Boltzmann  thermal  energy,  is  the  free  space  permeability,  and  Af 
is  the  average  domain  density.  The  saturation  magnetization  Ms  can  often  be  specified  from 
measured  data  or  a  priori  knowledge  about  the  material.  The  parameter  a,  on  the  other  hand, 
must  be  estimated  through  either  a  least  squares  fit  to  data  or  adaptive  parameter  estimation 
techniques  since  A f  is  unknown.  Due  to  its  dependence  on  the  temperature  T,  this  parameter 
may  vary  during  operation,  in  which  case,  adaptive  estimation  and  control  techniques  may  be 
preferable. 
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Hysteresis  losses  are  incorporated  through  the  quantification  of  energy  required  to  bend  and 
translate  domain  walls.  As  developed  in  [2],  an  irreversible  component  to  the  magnetization 
Mirr  can  be  computed  through  the  consideration  of  energy  dissipation  due  to  pinning  and 
unpinning  of  domain  walls  at  inclusions  in  the  material.  This  yields  the  expression 

dMirr  _  Man  Mj/Tr 

dH  ~  kd-a  (Man  -  Mirr )  1  ’ 

for  the  differential  susceptibility  where  k  =  .2m^\i-c)  ’  anc^  (p)  *s  average  density  of  pinning 
sites,  ( e j,.)  is  the  average  energy  for  180°  walls,  c  is  a  reversibility  coefficient,  and  m  is  the 
magnetic  moment  of  a  typical  domain.  Through  its  definition,  k  provides  a  measure  for  the 
average  energy  required  to  break  a  pinning  site.  The  parameter  5  is  defined  to  have  the  value 
+1  when  dH  >  0  and  —1  when  dH  <  0  to  guarantee  that  pinning  always  opposes  changes 
in  magnetization.  In  applications,  5  can  be  directly  determined  from  the  magnetic  field  data 
while  k  is  identified  for  the  specific  transducer  and  operating  conditions. 

We  note  that  in  certain  regimes  near  saturation,  the  expression  (3)  can  produce  nonphys¬ 
ical  results.  Specifically,  when  dH  is  reversed  near  saturation,  the  expression  (3)  can  produce 
negative  values  of  the  differential  susceptibility  which  are  not  observed  in  experiments.  Fol¬ 
lowing  the  strategy  detailed  in  [3],  we  consider  the  initial  changes  in  magnetization  to  be  due 
to  the  reversible  relaxation  of  bulged  domain  walls  when  dH  is  reversed  and  enforce  d^r  —  0. 
This  produces  the  more  physically  realistic  expression 


dMirr  _ ~  Man  I\IjrT 

dH  kd  —  a  ( Man  —  Mirr ) 


(4) 


where 


5  = 


1  ,  {dH  >  0  and  M  >  Man }  or  {dH  <  0  and  M  <  Man } 
0  ,  otherwise . 


The  reversible  magnetization  quantifies  the  degree  to  which  domain  walls  bulge  before 
attaining  the  energy  necessary  to  break  the  pinning  sites.  As  derived  in  [2],  to  first  approxi¬ 
mation,  the  reversible  magnetization  is  given  by 


Mrev  =  c(Man  -  Mirr) . 


(5) 


The  reversibility  coefficient  c  can  be  estimated  from  the  ratio  of  the  initial  and  anhysteretic 
differential  susceptibilities  [3]  or  through  a  least  squares  fit  to  data. 

The  total  magnetization  is  then  given  by 


M  —  Mrev  T  Mirr  (6) 

with  Mirr  and  Mrev  defined  by  (4)  and  (5)  and  the  anhysteretic  magnetization  given  by 
(1)  or  (2).  To  provide  a  relation  which  facilitates  inversion,  it  is  advantageous  to  express 
the  output  magnetization  as  a  function  of  the  input  field.  When  the  Ising  spin  model  is 
used  to  characterize  the  anhysteretic  magnetization,  the  expressions  (2),  (4)  and  (5)  can  be 
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consolidated  to  yield 


dM  dMirr  dMan 

m  =  ~  c)~dir + c-dH~ 

r(1_  s  Man  -  Mirr  dMan 

1  C)kS  -  a{Man  -  Mirr )  dH 

5(Man  -  M )  dMan 

k5  -  ^-c(Man  -  M)  °  dH 

_  l Ktaah  (IZZI  ~  M]  |  cM-,rr^(H  +  aM\  (l  I  ydM\ 

'  kS-  a  [M.tanh  (M!)  -  m]  +  o  l  a  ){ 1  +  a  dH  ) 


where  a  =  y^.  The  magnetization  at  a  given  field  level  is  then  specified  by  the  solution  to 
the  differential  equation 


dM 

~dH 


M ) 


(7) 


M(H0)  =  M0 


where 


F{H,  M )  = 


1  -  ^^sech2  f S±sM\ 

a  \  a  J 


5 

Mstanh  | 

(H+aM) 

\  a  J 

|  — 

M 

k5  — 

a 

Mstanh  | 

(  H+aM) 
{  a  J 

1  - 

M 

(8) 


cMs  2  (  H  +  aM 
- sech 


a 


a 


If  one  employs  the  Langevin  expression  (1)  rather  than  the  Ising  spin  relation  (2)  for  the 
anhysteretic,  the  function  T  is  given  by 


T{H,  M ) 


l 

r  5  [msc  (H+aM)  m] 

i  |  cMsa 
a 

(fflj-)  -  (5rihf)2' 

1  k5  a[Ms£(H+aM )  M] 

cMs 

a 


csclr 


H  +  aM 


H  +  aM 


(9) 


where  the  Langevin  function  is  defined  by 

C(z)  =  coth(z) - . 

z 

The  relation  (7)  with  T  given  by  (8)  or  (9)  is  then  employed  when  quantifying  the  actuator 
hysteresis  JV(-).  In  both  cases,  the  parameters  5  and  5  can  be  computed  directly  from  the 
measured  field  data  or  the  computed  values  of  the  magnetization  and  anhysteretic  magnetiza¬ 
tion.  The  saturation  magnetization  Ms  can  often  be  determined  directly  from  the  measured 
data  or  a  priori  knowledge  of  the  material  behavior.  Hence  only  the  parameters  a,  a,  c  and 
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k  must  be  determined  through  a  least  squares  fit  to  data  or  adaptive  parameter  estimation 
techniques. 

Finally,  we  note  that  the  time-dependent  dynamics  of  the  magnetization  can  be  specified 


through  the  chain  rule 


dM  _ .  TT  dH 
—  =  F(H,M)—. 
dt  '  dt- 


This  relation  should  be  employed  only  for  quasistatic  processes,  however,  since  this  model 
neglects  the  elastic  properties  of  the  material  as  well  as  time-dependent  loss  mechanisms. 


3  Inverse  Hysteresis  Model 

The  hysteresis  model  in  the  form  (7)  is  amenable  to  inversion  through  consideration  of  a 
complementary  differential  equation.  For  actuators  in  which  the  anhysteretic  is  modeled  by 
the  Ising  spin  model,  the  inverse  is  specified  through  solution  of  the  differential  equation 

dM-1  _  1 

dH  ~  T(M  l,H) 

M-\H0)  =M0-1 

with  T  defined  by  (8).  Similarly,  employment  of  (9)  provides  the  inverse  when  the  Langevin 
model  is  used  to  characterize  the  anhysteretic  magnetization. 

It  can  be  readily  verified  that  the  inverse  defined  through  (10)  satisfies 

M(M-l(H))  =  M-\M(H ))  =  H 

for  an  input  field  H.  In  the  control  nomenclature  of  Section  1,  the  definition  (10)  can  be 
used  to  provide  the  exact  inverse  fV_1(-)  if  the  parameters  Ms,  a,  a,  c  and  k  are  known  exactly. 
For  the  examples  in  this  paper,  it  is  assumed  that  suitable  parameters  have  been  obtained 
through  a  least  squares  fit  to  data  and  the  exact  inverse  is  employed.  Similarly  it  provides  an 
approximate  inverse  iV_1(-)  which  is  suitable  for  adaptive  parameter  estimation  and  control 
methods  if  the  parameters  are  unknown  or  slowly  change  due  to  changing  operating  conditions. 

Example  1. 

In  control  applications,  one  typically  employs  the  computed  inverse  as  a  filter  to  com¬ 
pensate  for  the  physical  hysteresis  and  nonlinearities  produced  by  the  actuator.  Hence  one 
considers  the  operation 

u  =  N(N-1(up )) 

where  up  is  the  specified  control  in  the  absence  of  actuator  hysteresis  and  nonlinearities.  This 
example  illustrates  the  inversion  process  with  the  hysteresis  specified  by  N  =  M. 

The  model  (7)  and  inverse  (10)  were  constructed  with  the  parameters  Ms  =  0.280  C/m2, 
a  =  4.2  x  105C/m2,k  =  3.7  x  10 5C/m21a  =  1.5  x  106,c  =  0.65.  A  1  Hz  piecewise  linear 
signal  up  having  a  maximum  magnitude  of  0.25  V/m  was  employed  as  input.  Such  quasistatic 
signals  are  commonly  employed  in  material  characterization  experiments.  A  first-order  forward 


(10) 
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difference  method  (Euler’s  method)  was  used  to  numerically  integrate  the  differential  equations 
(7)  and  (10).  Due  to  the  high  sample  rates  in  typical  experiments,  the  method  provides 
sufficient  accuracy,  efficiency  and  stability  for  both  numerical  and  real  time  experimental 
implementation. 

The  input  signal  uP  and  output  signal  u  =  M{M^l{up))  are  plotted  in  Figure  3  while  the 
inverse  function  v  =  M~l(u,p)  is  plotted  in  Figure  4.  A  comparison  of  up  and  M(M^1(up)) 
illustrates  the  qualitative  accuracy  of  the  inversion  process  while  numerical  computations 
indicate  that  the  absolute  errors  introduced  in  the  process  have  magnitude  less  than  10~14. 
While  such  accuracy  obviously  cannot  be  expected  when  the  finite  dimensional  inverse  iV_1(-) 
is  employed  in  an  infinite  dimensional  process  JV(-),  the  maintenance  of  these  accuracy  levels 
under  approximation  does  illustrate  the  stability  of  the  method. 


Figure  3.  (a)  Signal  up(t)  input  to  the  inverse  compensator;  (b)  Signal  u(t)  =  M(M  1{up(t))) 
applied  to  the  plant. 


Figure  4.  Inverse  function  M  i(up). 
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4  Structural  Application  and  Control  Example 


To  illustrate  the  incorporation  of  the  actuator  hysteresis  model  in  a  structural  model  and 
indicate  the  utility  of  the  inverse  compensator  in  linear  control  design,  we  consider  a  can¬ 
tilevered  thin  beam  with  end-mounted  magnetostrictive  actuators  as  depicted  in  Figure  5. 
Spatially  uniform  forces  drive  the  beam  while  diametrically  out-of-phase  currents  to  the  actu¬ 
ators  generate  bending  moments  which  attenuate  transverse  beam  vibrations.  This  setup  has 
been  employed  in  experiments  designed  to  ascertain  capabilities  and  properties  of  Terfenol-D 
transducers  [6]  and  it  provides  a  prototype  for  illustrating  the  attributes  of  the  inverse  com¬ 
pensator.  For  further  details  concerning  the  construction  of  the  actuator,  the  reader  is  referred 
to  [6]. 

For  this  example,  we  will  employ  an  LQR  feedback  control  law.  However,  the  inverse 
compensator  is  equally  effective  for  linear  feedforward  or  adaptive  control  methods  as  dictated 
by  the  application. 

For  modeling  purposes,  the  beam  is  assumed  to  have  length  £,  width  b,  and  thickness  h. 
The  density,  Young’s  modulus,  Kelvin- Voigt  damping  coefficient  and  air  damping  coefficient 
for  the  beam  are  denoted  by  pb,  Eb,  cub  and  7,  respectively.  The  cross-sectional  area  of  the 
Terfenol  rod  is  denoted  by  Amag  while  the  Young’s  modulus  and  damping  coefficient  for  the 
Terfenol  rod  are  denoted  by  EH  and  c^.  The  length  and  width  of  the  connecting  bar  are 
denoted  by  £r  and  br.  respectively,  while  the  bar  density  is  given  by  pr.  Finally,  the  transverse 
beam  displacement  is  given  by  w  while  g(t.  x )  denotes  an  exogenous  surface  force  to  the  beam. 

Moment  and  force  balancing  yields  the  strong  form  of  the  Euler-Bernoulli  equations 


0  <  x  <  £ 
t  >  0 


w(t10)  =  ^(t10)  =  0  ^ 

ux  t  >  0 

Mint(t,£)  =  d^nt  (M)  =  0 

along  with  appropriate  initial  conditions,  as  a  model  for  characterizing  the  transverse  beam 
dynamics.  As  detailed  in  [7],  the  composite  density  and  internal  bending  moment  are  given 
by 

p(x)  =  pbhb  +  2prbr£rXrod{x ) 

d^w  d^w 

Mint(t,x)  =  EI(x)^(t,x)  +  cDI-Q^(t,x) 
where  the  characteristic  function  Xrod  delineates  the  location  of  the  rods  and 

E  h?b 

EI 0)  =  — +  2 AmagEH  (h/2  +  £rf  Xrod(x) 

cDI(x )  =  CDb^  b  +  2AmagCp  (h/2  +  £rf  Xrod(x) . 

To  specify  the  external  moments  A4mag,  it  is  necessary  to  quantify  the  strains  generated 
by  an  input  field  H  to  the  actuator.  While  the  strains  depend  upon  a  variety  of  material  prop- 
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Figure  5.  Cantilever  beam  with  magnetostrictive  actuators.  Uniform  force  inputs  are  de¬ 
picted  above  the  beam  while  the  measurement  point  is  indicated  by  the  lower  arrow. 


erties  including  the  crystalline  anisotropy,  to  first  approximation  they  can  be  specified  by  the 
bulk  magnetostriction 

9  \ 

m  =  (12) 

where  Xs  denotes  the  saturation  magnetostriction  and  M  is  the  magnetization  given  by  (7). 
The  combination  of  (7)  and  (12)  provides  a  relation  between  the  input  field  H  and  positive 
strains  generated  by  the  actuator.  To  obtain  bi-directional  strains,  the  weighted  magnetization 
2 M(t)Ms  is  used  to  bias  to  the  center  of  the  magnetic  range.  This  yields  the  strains 

9  \ 

e(t)  =  2ii^[M2(t)  +  2M(t)Ms] 


along  with  the  expression 

Mmag(t,x)  =  ICM[M2(t )  +  2  M(t)Ms]Xrod(x) 


for  the  moment  produced  by  diametrically  out-of-phase  Terfenol-D  rods.  Here  the  constant 
KM  is  given  by  (3A s/M2)AmagEH  (h/2  +  £r)2 . 

To  obtain  a  weak  form  of  the  model,  we  take  the  state  to  be  the  displacement  w  in  the 
state  space  X  =  L2(0,£)  with  the  inner  product 


ri 

(<M)x  =  /  pHdx. 

Jo 


The  space  of  test  functions  is  taken  to  be  V  =  H'l(0,£ )  =  {(f)  e  H2(0,£)  \  (f)(0)  =  <^'(0)  =  0} 
with  the  inner  product 

{(f),ip)v  =  f  E  1(f)" 'll;"  dx  . 

Jo 

A  weak  form  of  the  model  is  then  given  by 

/  pw(f)dx+  /  ryw(f)dx+  /  M.intcf)"  dx  =  /  Aimag4>"  dx  +  /  g(f>dx  (13) 

Jo  Jo  Jo  Jo  Jo 

for  all  (j)  €  V.  It  is  in  this  form  that  we  develop  the  approximation  method  and  formulate  the 
control  problem. 
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To  approximate  the  dynamics  of  (13)  in  a  manner  amenable  to  control  formulation,  we 
employ  a  Galerkin  discretization  in  the  spatial  variable  to  obtain  a  semidiscrete  ODE  system 


■th 


in  time.  Specifically,  the  spatial  basis  is  taken  to  be  where  4>j{x)  denotes  the  j 

cubic  5-spline  modified  to  satisfy  the  fixed  left  boundary  condition.  Approximate  solutions 

(14) 


m+ 1 

wm(t,x)  =  wAMjip) 
j  1 


are  then  considered  in  the  subpace  Vm  =  span{^j}.  To  obtain  a  vector  ODE  system,  the 
infinite  dimensional  system  (13)  is  restricted  to  Vm  and  posed  in  first-order  form  to  yield 


y(t)  =  Ay(t)  +  B(u(t ))  +  G(t) 
2/(0)  =  2/o  • 

The  component  system  matrices  have  the  form 

0  / 

Q-lK  Q-'C 


A 


B(u(t ))  =  [ M2(u{t ))  +  2  M(u(t))Ms] 


0 

Q  'B 


0 

Q^g(t) 

where  y(t)  =  [w^t),  ■  ■  ■ ,  wm+1(t),  w^t),  ■  ■  ■ ,  wm+1(t)]  and 

K,m 


[Q\ij 

[K]ij 


f  p<pi4>j  dx  [B\i  =  K,m  f  (()"  dx 

JO  J  mag 

J  EI(t>"<t>”  dx  [g(t)]i  =  g{t ,  x)<f>i  dx 


[C]ij  =  J  cdI4>"4>"  dx  . 


Note  that  u(t)  =  H(t )  denotes  the  control  input  to  the  system. 

From  (16),  we  note  that  the  hysteresis  input  operator  has  the  form 

N(H(t))  =  B(H(t))  =  [ M2(H(t ))  +  2 M(H(t))Ms]  b 

where  b  denotes  the  2 (m+  1)  vector 

t _ r  o 

1  Q  lB 

The  corresponding  inverse  is 

N-l(H(t))  =  M-1  (-M,  +  v/:V/(  -  H(l)  ) 


(15) 


(16) 


(17) 


(18) 
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where  AW1  denotes  the  solution  to  the  differential  equation  (10). 

The  specific  form  of  the  optimal  control  is  dependent  on  the  nature  of  the  exogenous  force 
g  in  (11)  or  (13)  and  resulting  disturbance  term  G  in  (15).  For  the  discussion  which  follows, 
we  consider  exogenous  forces  which  are  periodic  with  fundamental  period  r.  Such  forces  arise 
when  modeling  systems  driven  by  oscillating  or  rotating  components. 

We  now  consider  two  techniques  for  formulating  a  linear  optimal  control  problem.  In  the 
first  case,  we  simply  linearize  the  operator  B(H )  about  a  given  field  level.  This  technique 
is  commonly  employed  when  magnetostrictive  transducers  are  employed  at  low  drive  levels 
with  magnetic  biases  and  is  applicable  only  when  the  transducer  is  operated  in  the  nearly 
linear  range.  The  second  technique  is  based  upon  the  filtering  of  control  inputs  using  the 
inverse  compensator.  This  control  method  is  applicable  throughout  the  operational  range  of 
the  transducer. 

To  specify  the  control  input  u(t)  =  H (t)  in  the  first  case,  we  will  linearize  about  the 
coercivity  value  uc  at  which  M(uc )  =  0.  In  this  case,  the  approximate  linear  operator  B  is 


B  =  2Ms^-{uc)b 
du 


(19) 


where 


dM 

du 


5{1  -  c) 


Man  ~  Mirr 
kS  -  a[Man  -  Mh 


+ 


cAL 


csclT 


H, 


eff 


a 


+ 


for  the  anhysteretic  model  (1). 

We  then  consider  the  minimization  of 


«/(«) =  \L  [yT(t')Qy(t') +  uT(t)Ru(t )] dt 


H, 


eff , 


(20) 


subject  to 

y(t)  =  Ay(t)  +  Bu(t)  +  G(t) 

3/(0)  =  y(r)  • 

Here  Q  and  R  respectively  denote  the  matrices  used  to  weight  the  state  and  control  input. 
Under  typical  stabilizability  and  detectability  assumptions  (e.g.,  see  [8]),  the  optimal  control 
for  (20)  is  given  by 

u(t)  =  -R^1BT[Uy(t)  -  r(t )]  (22) 


where  n  solves  the  algebraic  Riccati  equation  HTn  +  n.4  —  IiBRrlBTB  +  Q  =  0 .  The  per¬ 
turbation  variable  r(t)  e  ]R2^m+1^  is  obtained  through  solution  of  the  periodic  perturbation 
system 


r(t)  =  -[A  -  BRr1BTU]Tr(t)  =  UG(t) 
r(0)  =  r(r) . 


(23) 


Note  that  in  applications,  (23)  is  typically  integrated  backward  in  time  from  r  or  forward  in 
time  from  —  r. 

As  will  be  demonstrated  in  subsequent  examples,  the  application  of  the  control  u(t)  spec¬ 
ified  by  (22)  into  the  original  system  (15)  is  effective  at  low  input  levels  but  fails  drastically 
at  high  drive  levels  as  a  result  of  unincorporated  time  delays  due  to  the  hysteresis. 
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A  control  based  on  inverse  compensation  is  obtained  by  employing  the  input  operator 

B  =  b 


where  b  is  given  by  (18).  This  yields  the  resulting  input  up.  To  obtain  the  optimal  control  for 
the  original  system  (15),  the  input  up  is  filtered  through  the  inverse  to  obtain 

v(t)  =  M  1  (m,  +  y/M*  +  up{t)  )  .  (24) 

We  note  that  while  this  requires  the  online  integration  of  (10)  to  obtain  M_1,  the  operation  is 
significantly  less  expensive  than  the  online  integration  of  (23)  to  obtain  r(t )  or  the  integration 
required  to  obtain  state  estimates  yc(t )  when  full  state  information  is  unavailable. 


Example  2. 

To  illustrate  the  performance  and  limitations  of  the  two  techniques,  we  consider  a  beam 
of  length  l  =  ,4573  m,  width  b  =  ,0203  m  and  thickness  h  =  ,0016  m  with  the  material 
parameters  for  the  beam  and  actuators  summarized  in  Table  1.  It  should  be  noted  that  the 
beam  parameters  are  consistent  with  typical  values  for  aluminum  laboratory  beams  while 
the  Terfenol  parameters  are  within  the  range  obtained  for  model  fits  to  an  experimental 
transducer  [4],  For  this  choice  of  beam  parameters,  the  first  two  natural  frequencies  for  the 
system  occur  at  6.1  Hz  and  38.3  Hz.  To  account  for  the  effects  of  parameter  discontinuities 
due  the  actuators  and  damping  in  the  system,  it  was  necessary  to  obtain  these  values  through 
a  fast  Fourier  transform  (FFT)  of  time  domain  data  resulting  from  a  simulated  impact  to  the 
beam  (it  is  not  possible  to  obtain  analytic  expressions  through  separation  of  variables).  The 
driving  frequency  in  the  numerical  examples  was  chosen  close  to  but  not  exactly  concurrent 
with  these  natural  frequencies.  For  this  example,  the  choice  m  =  12  was  sufficient  for  resolving 
the  beam  dynamics  in  the  frequency  range  considered  and  the  reported  results  were  obtained 
with  m  =  16. 

A  spatially  uniform  exogenous  force 

g(t,  x )  =  g0  [sin(147rf)  +  sin(267rf)] 

was  applied  throughout  the  time  interval  [0,  2.5]  to  simulate  a  periodic  pressure  field  on  the 
beam.  The  uncontrolled  trajectories  at  the  point  x  =  3^/5  for  go  =  l  are  plotted  in  Figures  6b 
and  7a  while  the  trajectories  obtained  with  g0  =  100  are  plotted  in  Figures  6d  and  7b.  Both 
cases  exhibit  a  beat  phenomenon  due  to  the  close  proximity  of  the  7  Hz  driving  frequency 
with  the  6.1  Hz  natural  frequency  of  the  beam. 

We  first  consider  the  performance  of  the  control  based  on  the  linearized  operator  B  defined 
in  (19).  The  trajectories  which  result  when  the  input  u  given  by  (22)  is  fed  back  into  the 
hysteretic  system  (15)  at  time  t  =  0.45  are  plotted  in  Figures  6b  and  6d.  The  relationship 
between  the  input  magnetic  field  H  =  u  and  the  resulting  magnetization  M  for  the  two  cases 
are  given  in  Figures  6a  and  6c.  The  results  for  g0  =  1  illustrate  that  the  method  is  very  effective 
at  low  drive  levels  where  the  linear  model  is  accurate.  At  the  high  drive  levels  at  which  the 
actuators  are  advantageous,  however,  the  input  B(u )  introduces  significant  hysteresis  which 
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acts  as  a  time  delay  to  the  system.  The  result  is  a  loss  in  control  authority  which  is  sufficiently 
severe  to  produce  controlled  beam  trajectories  which  are  larger  in  magnitude  at  certain  times 
than  the  uncontrolled  trajectories. 

The  failure  of  the  control  based  on  the  linearized  system,  when  employed  at  high  drive 
levels,  motivates  the  consideration  of  a  control  based  on  inverse  compensation.  In  this  case,  the 
operator  B  —  b  yields  an  input  up  which  is  filtered  by  the  inverse  compensator  (24)  before  input 
to  the  system.  As  illustrated  in  Section  3,  the  resulting  input  u  to  the  plant  is  then  precisely 
the  control  up  specified  in  the  absence  of  hysteresis  and  nonlinear  dynamics.  The  controlled 
trajectories  for  g0  =  1  and  g0  =  100  are  respectively  plotted  in  Figures  7a  and  7b.  In  this  case, 
there  is  no  loss  of  control  authority  at  high  drive  levels,  thus  illustrating  that  this  approach 
provides  a  significant  advantage  over  the  control  design  based  on  the  linearized  system.  Finally, 
a  comparison  with  the  examples  in  [9]  illustrates  that  the  performance  of  the  controller  based 
on  inverse  compensation  is  comparable  to  that  obtained  when  the  functional  (20)  is  minimized 
subject  to  the  original  system  (15)  but  at  a  fraction  of  the  computational  cost. 


(c)  (d) 


Figure  6.  Feedback  of  the  linear  law  (22)  into  the  nonlinear  system  (15).  Relationship 
between  magnetic  field  H  and  magnetization  M;  (a)  go  =  1  and  (c)  g0  =  100.  Uncontrolled 
( - )  and  controlled  (  )  beam  trajectories  at  x  =  3£/5]  (b)  go  =  l  and  (d)  g0  =  100. 
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(a)  (b) 

Figure  7.  Uncontrolled  beam  trajectory  ( - )  at  the  point  x  =  3^/5  and  controlled  trajectory 

(  )  obtained  with  the  filtered  input  (24);  (a)  g0  =  1  and  (b)  g0  =  100. 


Beam 

Actuator 

Terfenol 

Eb  =  7.0861  x  1010  N/m2 

EH  =  7.0  x  1010  N/m 2 

a  =  7105  A/m 

pb  =  2863  kg/m 3 

pr  =  8524  kg/m 3 

k  =  7002  A/m 

cDb  =  9.3663  x  105  Ns/m2 

tr  =  .0254  m 

a  =  .007781 

7  =  .013  Ns/nri2 

br  =  .002  m 

Amag  =  .0064  m2 

c  =  0.3931 

Ms  =  1.3236  x  105  A/m 
Xs  =  9.96  x  10-4 

Table  1.  Parameters  for  the  beam  and  Terfenol  transducer. 


5  Concluding  Remarks 

This  paper  addresses  the  development  of  a  modeling  strategy  and  corresponding  model-based 
inverse  compensator  for  a  class  of  ferromagnetic  transducers  including  magnetostrictive  ac¬ 
tuators.  Control  designs  based  on  linearized  system  models  are  typically  ineffective  in  such 
regimes  since  they  do  not  provide  the  capability  for  incorporating  hysteresis  and  subsequent 
time  delays.  The  result  is  a  severe  loss  in  control  authority  and  potential  destabilization  of 
the  system.  To  address  this,  we  consider  a  model  inverse  which  provides  a  transform  so  that 
linear  inputs  produce  linear  outputs  when  applied  to  the  actuator. 

The  models  are  based  on  the  Jiles-Atherton  theory  which  quantifies  the  energy  required 
to  bend  and  translate  domain  walls  pinned  at  inclusions  in  the  materials.  The  construction 
of  the  models  can  be  summarized  in  the  following  four  steps:  (1)  Determine  the  effective 
field  at  a  point  in  the  material,  (2)  Quantify  the  anhysteretic  magnetization  through  Boltz¬ 
mann  principles,  (3)  Quantify  the  irreversible  magnetization  due  to  domain  wall  translations 
and  (4)  Characterize  the  reversible  magnetization  due  to  domain  wall  bending.  The  total 
magnetization  is  then  the  sum  of  the  irreversible  and  reversible  components. 
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The  consolidated  model  quantifying  the  magnetization  in  ferromagnetic  materials  is  in  the 
form  of  an  initial  value  problem  which  can  be  marched  in  time  to  characterize  the  hysteresis 
loop.  A  complementary  differential  equation  then  provides  the  inverse  transforms  for  the 
actuators.  The  performance  of  this  compensator  and  its  employment  in  LQR  control  design  are 
illustrated  through  numerical  examples.  It  is  noted  that  while  the  determination  of  the  inverse 
requires  the  numerical  approximation  of  a  differential  equation,  the  operation  is  significantly 
less  expensive  than  the  online  integration  required  to  compute  perturbation  variables  or  update 
state  estimates.  Hence  the  inverse  compensator  appears  feasible  for  real-time  experimental 
implementation. 

Finally,  we  point  out  that  while  not  discussed  here,  similar  models  have  been  developed 
to  quantify  hysteresis  in  certain  ferroelectric  materials  [10].  Using  techniques  similar  to  those 
described  here,  the  use  of  these  models  provides  the  capability  for  compensating  for  hysteresis 
observed  in  relaxor  ferroelectric  materials  at  low  temperatures  and  piezoceramic  actuators  at 
high  drive  levels.  The  comprehensive  development  of  inverse  compensators  based  on  ferro¬ 
electric  domain  wall  theory  is  under  investigation  and  will  be  reported  in  a  future  work. 
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